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DYNAMIC  FRACTURE  IN  TEMPERATURE  SENSITIVE  MATERIALS 
WITH  COHESIVE  ZONES:  A  DISCONTINUOUS  GALERKIN  APPROACH 

GRANT  NO.  F496 20-02- 1-0318 

Francesco  Costanzo 

Engineering  Science  and  Mechanics  Department 
212  Earth  and  Engineering  Sciences  Building 
The  Pennsylvania  State  University 
University  Park,  PA  16802 

Phone:  (814)  863-2030;  E-Mail:  costanzo@engr.psu.edu 


1.  Objectives 

This  project  intends  tp  expand  the  understanding  of  the  role  of  temperature  in  controlling  the 
dynamic  failure  behavior  of  advanced  materials  subject  to  combined  thermo-mechanical  loading. 
These  objectives  will  be  achieved  by  improving  the  continuum-based  modeling  of  the  fracture 
properties  of  materials  in  conjunction  with  the  formulation  and  development  of  a  corresponding 
numerical  approach  for  the  solution  of  the  resulting  governing  equations. 

The  goal  of  the  proposed  modeling  effort  is  the  inclusion  of  temperature  dependence  in  the 
description  of  the  fracture  behavior  of  the  material.  In  addition  to  modeling  the  bulk  material  as  a 
fully  coupled  linear  thermo-elastic  solid,  this  goal  will  be  achieved  by  using  rate  and  temperature 
dependent  cohesive  zone  (CZ)  models  to  account  for  the  highly  localized,  nonlinear,  and  softening 
behavior  of  the  material  ahead  of  the  propagating  crack.  In  this  project,  the  cohesive  stresses 
will  be  assumed  to  be  physically-based  constitutive  functions  of  the  opening  displacement,  opening 
displacement  rate  and  temperature.  Finally,  the  CZ  will  be  assumed  to  fail  in  two  basic  ways:  (?) 
by  achieving  a  critical  crack  opening  displacement;  and/or  («)  by  experiencing  a  critical  value  of 
cohesive  stress.  Both  the  critical  crack  opening  displacement  and  the  critical  cohesive  stress  will 
be  assumed  to  be  functions  of  temperature. 

The  primary  objective  of  the  proposed  numerical  development  is  the  formulation  of  a  high 
accuracy  and  unconditionally  stable  solution  scheme  for  the  combined  parabolic/hyperbolic  prob¬ 
lem  describing  dynamic  fracture  in  a  linear  thermo-elastic  solid.  Adaptivity  will  be  a  primary 
feature  of  the  proposed  numerical  scheme  as  it  will  be  is  crucial  for  the  accurate  quantification 
of  the  microscopic  features  that  are  nucleated  during  propagation  as  these  features  are  a  means 
of  comparison  between  theory  and  experiments.  Said  numerical  approach  will  be  based  on  the 
discontinuous  Galerkin  finite  element  method  (DG  FEM).  The  DG  FEM  has  been  shown  to  be 
extremely  effective  in  the  solution  of  both  parabolic  and  hyperbolic  problems  in  the  presence  of 
moving  discontinuities. 

To  provide  an  alternative  strategy  for  the  assessment  of  the  fracture  properties  resulting  from  a 
given  choice  of  the  CZ  constitutive  law,  and  to  provide  a  tool  to  test  the  efficacy  of  the  proposed  DG 
FEM  development,  the  proposed  research  also  includes  the  use  of  a  different  numerical  technique 
previously  developed  by  the  proposer  and  a  collaborator.  This  technique,  based  on  the  notions 
of  Neumann-to-Dirichlet  map  and  product  integration,  is  highly  accurate  but  only  applicable  to 
problems  with  idealized  geometry  (e.g.,  dynamic  propagation  of  a  semi-infinite  crack  along  a  planar 
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bi-material  interface).  The  aforementioned  technique  will  be  expanded  and  refined  for  the  purpose 
of  creating  an  effective  benchmarking  tool  for  the  proposed  DG  FEM. 


2.  Status  of  Effort  and  Accomplishments 

This  project  begub  July  1,  2002  and  will  expire  on  November  30,  2004. 

During  the  last  two  years,  the  primary  focus  of  the  this  effort  has  been  the  creation  of  a 

reliable  computational  capability  for  the  solution  of  elasto-dynamic  moving  boundary  problems. 

Specifically,  the  PI,  along  with  a  student, 

1.  has  formulated  a  space-time  discontinuous  Galerkin  finite  element  (DGFEM)  formulation  that 
has  been  shown  to  be  unconditionally  stable  and  capable  of  adaptive  mesh  refinement; 

2.  applied  this  formulation  to  the  solution  of  various  model  problems  in  isothermal  elasto- 
dynamics  as  well  as  isothermal  elasto-dynamic  phase  transition; 

3.  begun  the  extension  of  the  DGFEM  formulation  to  fully  coupled  thermo-elastic  problems  along 
with  preliminary  numerical  results;  and 

4.  begun  the  creation  of  a  numerical  procedure  for  the  solution  of  double  integral  equations 
deriving  from  the  reformulation  of  fully  transient  Mode  III  dynamic  fracture  problems  in 
thermo  elastic  media  with  cohesive  zones. 

2.1.  DGFEM:  early  developments 

The  choice  of  a  discontinuous  Galerkin  method  for  the  solution  of  dynamic  fracture  problems  was 

motivated  by  several  reasons: 

(a)  The  fragmentation  processes  occurring  ahead  of  a  running  crack  tip  are  intrinsically  discon¬ 
tinuous  in  both  space  and  time.  However,  second  order  accurate  numerical  methods  with  no 
artificial  viscosity  are  likely  to  predict  severe  high  frequency  spurious  oscillations  when  the 
boundary  data  are  non-smooth.  Clearly,  one  could  resort  to  methods  which  introduce  artificial 
viscosities  (see,  for  example,  the  discussion  on  the  /?- Newmark’s  method  by  Hughes,  1987  or 
Johnson,  1987).  However,  this  practice  is  questionable  in  the  present  context  due  to  the  fact 
that  the  CZ  constitutive  equations  are  rate  dependent  and  it  would  be  extremely  complicated 
to  distinguish  the  rate  effects  due  to  the  physical  behavior  of  the  near  crack  tip  environment 
from  the  behavior  induced  by  the  choice  of  the  numerical  solution  scheme.  In  addition,  one 
must  remember  that  high  frequency  crack  tip  velocity  oscillations  are  one  of  the  physical  phe¬ 
nomena  that  one  wishes  to  be  able  to  predict.  The  PI  intended  to  formulate  a  solution  scheme 
free  from  the  suspicion  that  high  frequency  oscillations  are  induced  by  the  solution  strategy 
itself.  Furthermore,  it  should  be  remembered  that  no  rigorous  studies  have  been  done  which 
discuss  error  estimation  in  dynamic  fracture  problems  and  therefore  it  is  very  difficult  to  as¬ 
sess  the  stability/accuracy  properties  of  current  solution  schemes.  By  contrast,  as  shown  by 
Hughes  and  Hulbert  (1988),  Hulbert  and  Hughes  (1990),  and  Johnson  (1993)  (see  also  French, 
1993;  Hulbert,  1992a, c)  second  order  hyperbolic  problems  can  be  given  a  DG  formulation  which 
leads  to  unconditionally  stable  and  high  accuracy  algorithms  without  the  addition  of  artificial 
viscosities.  In  addition,  the  space-time  formulation  of  the  problem  makes  it  possible  for  one  to 
conduct  a  very  transparent  error  analysis. 
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(b)  A  DGFEM  in  dynamic  fracture  requires  that  the  meshing  of  the  space-time  domain  be  dynamic, 
that  is,  that,  the  FEM  grid  be  generated  along  with  the  rest  of  the  computation.  This  is  crucial 
if  one  wants  to  avoid  the  pre-meshing  (that  is,  the  a  priori  assigning)  of  the  crack  path  as  it  is 
done  in  any  gradual  node  release  algorithm,  the  latter  being  by  far  the  most  used  computational 
strategy  in  dynamic  fracture  to  date. 

(c)  The  third  important  reason  for  proposing  the  use  of  a  DGFEM  approach  is  that  the  present 
proposal  intends  to  explore  the  fully-coupled  thermo-mechanical  dynamic  fracture  problem. 
This  problem  is  notoriously  difficult  to  confront  numerically  due  to  the  fact  that  semi-discrete 
FEM  formulations  (i.e.,  FEM  schemes  which  discretize  the  space  domain  only)  of  the  heat 
equation  lead  to  stiff  initial  value  problems  for  the  time  component  of  the  problem.  In  the 
context  of  parabolic  problems  (such  as  the  heat  equation)  the  advantage  of  the  use  of  DGFEM 
has  long  been  demonstrated.*  Therefore  the  use  of  a  DGFEM  approach  to  the  problem  of 
interest  herein  seems  rather  “natural”  since  one  will  be  able  to  construct  solutions  schemes  for 
which  error  estimates  can  be  assessed  in  a  rational  and  unified  way  for  both  the  mechanical 

and  thermal  parts  of  the  problem,  thus  leading  to  the  possibility  of  adaptive  and  parallelizable 
solution  schemes. 

As  part  of  the  work  done  in  the  last  2  years  (although  some  of  the  work  was  begun  before  August 
1,  2002,  as  part  of  the  Ph.D.  thesis  of  Hao  Huang,  one  of  the  Pi’s  students),  the  PI  revisited  the 
DGFEM  formulation  proposed  Hughes  and  Hulbert  (1988);  Hulbert  and  Hughes  (1990)  and  adapted 
it  to  the  study  of  moving  boundary  problems  in  isothermal  elasto-dynamics.  This  extension  has 
been  reported  in  Huang  and  Costanzo  (2002),  in  which  basic  properties  of  the  extended  formulation 
were  discussed  and  along  with  a  number  of  results  concerning  a  ID  model  problem  with  dynamic 
solid/solid  phase  transition  as  depicted  in  Fig.  1,.  The  motion  considered  is  that  of  a  linear  elastic 


Figure  1:  Bi-material  bar  in  tension  fixed  at  both  ends. 

ID  bi-material  bar  subjected  to  an  initial  stretch.  The  two  halves  of  the  bar  are  given  the  same 
density  p0  but  different  stiffiiess,  k\  and  k2,  respectively,  with  ki  <  k2,  and  correspondingly  different 
wave  speeds,  ci  and  c2,  with  ci  <c2.  At  time  t  —  0  the  interface  separating  the  two  halves  of  the  bar 
is  made  to  move  in  a  prescribed  manner  such  that  the  less  stiff  part  of  the  bar  grows  at  the  expense 
of  the  stiffer  part.  Said  motion  is  therefore  analogical  to  a  damage  or  fracture  process  since  it  causes 
a  weakening  of  the  overall  structure.  The  results  reported  below  are  presented  in  non-dimensional 
form.  In  particular,  the  quantities  reported  in  the  graphs  have  been  non-dimensionalized  as  follows: 


x*  =  -  -  £2,  ,  _  1 

V  *  “  L  '  ~F2  ’ 


(1) 


There  is  vast  literature  on  the  use  of  discontinuous  Galerkin  FEM  in  the  solution  of  linear  and  nonlinear  parabolic 
problems  and  its  review  is  outside  the  scope  of  this  proposal.  Here  are  a  few  references  which  focus  on  the  heat 
equation:  Cockburn  and  Shu  (1998);  Eriksson  and  Johnson  (1991,  1995);  Eriksson  et  al.  (1985);  French  (1998  1999V 
(1996eT200of2b  ’  J°hnSOn  (198?)i  Makridakis  and  Babuska  (1997);  Schotzau  and  Schwab  (2000);  Shaw  and  Whiteman 
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Figure  3:  Distribution  of  the  axial  internal  force  at  various  time  instants  for  the  case  of  the  interface 
traveling  at  a  constant  acceleration  a.  Here  (aL)/c$  =  0.2,  c\/c2  —  0.5,  k\/k2  =  0.25,  and 
n*(t  —  0)  =  0.1.  Every  plot  has  the  same  horizontal  and  vertical  axes.  The  horizontal  axis 
represents  the  non-dimensional  position  x*  along  the  bar,  while  the  vertical  axis  represents  the 
non-dimensional  axial  force  n*.  The  space  time  domain  was  discretized  using  400  bi-quadratic 
isoparametric  space-time  elements. 
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where  x  and  x*  denote  the  position  and  the  non-dimensional  position  along  the  bar,  respectively,  t 
and  t*  denote  time  and  non-dimensional  time,  respectively,  n  and  n*  denote  the  bar’s  internal  axial 
force  and  its  non-dimensional  counterpart,  respectively.  The  first  set  of  results  concerns  the  case 
of  an  interface  traveling  along  a  prescribed  trajectory.  Two  cases  are  reported  here:  at  a  constant 
speed  case  and  a  constant  acceleration  case  reported  in  Fig.  2  and  Fig.  3,  respectively.  For  the 
constant  velocity  case,  it  should  be  noted  that  the  growth  process  modeled  is  discontinuous  in  two 
distinct  ways.  Firstly,  the  process  is  discontinuous  because  the  motion  of  the  interface  is  such  that 
the  interface  simply  jumps  from  rest  a  steady  state  regime.  Secondly,  the  process  is  discontinuous 
in  that  the  internal  axial  force  in  the  bar  suffers  jumps  across  the  propagating  acoustic  wave  fronts. 
Despite  these  discontinuity  sources  the  calculated  results  are  remarkably  close  to  those  given  by  the 
exact  solution,  especially  in  view  of  the  fact  that  the  quantity  n*  is  obtained  from  the  derivative  of 
the  primary  unknown,  the  latter  being  the  axial  displacement  of  the  bar. 


Figure  2:  Distribution  of  the  axial  internal  force  (stress)  at  various  time  instants  for  the  case 
of  the  interface  traveling  at  a  constant  speed  with  v/c2  —  0.1.  C1/C2  =  0.5,  fci/&2  =  0.25, 
and  n*(t  ~  0)  =  0.1.  Every  plot  has  the  same  horizontal  and  vertical  axes.  The  horizontal 
axis  represents  the  non-dimensional  position  x*  along  the  bar,  while  the  vertical  axis  represents 
the  non-dimensional  axial  force  n*.  Each  space-time  slab  was  discretized  using  400  bi-quadratic 
isoparametric  space-time  elements. 

To  quantify  the  accuracy  of  the  method,  convergence  rates  for  the  two  cases  reported  herein 
have  been  explicitly  computed  along  with  the  convergence  rate  for  a  problem  yielding  a  smooth 
solution.  In  particular,  as  suggested  in  Hulbert  and  Hughes  (1990),  the  solution  for  the  case  of  a 
homogeneous  bar  fixed  at  both  ends  with  a  quiescent  past  is  subjected  to  a  given  initial  displacement 
field  coinciding  with  the  displacement  distribution  given  by  the  first  harmonic  of  the  system.  For 
this  test  case  the  convergence  rate  is  3,  as  theoretically  predicted  by  Hulbert  and  Hughes  (1990). 
Clearly,  the  convergence  rate  is  negatively  affected  by  the  loss  of  smoothness  in  the  solution  of 
the  moving  interface  problem.  However,  as  Fig.  4  indicates,  the  convergence  rates  for  the  moving 
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interface  problems  discussed  herein  are  still  quite  satisfactory. 


Convergence  Rates 


log(l//i) 


Figure  4:  Convergence  rates  for  a  problem  with  smooth  solution  along  with  those  of  the  solid-solid 
phase  transition  problems  preciously  discussed. 

The  proposed  DGFEM  has  also  been  tested  by  the  PI  in  the  solution  of  a  few  fracture  model 
problems  with  excellent  results,  which  have  been  reported  in  Huang  and  Costanzo  (2004).  Problems 
with  sharp  cracks  as  well  as  with  a  Dugdale  zones  have  been  considered.  For  the  purpose  of 
comparison,  we  report  the  energy  release  rate  computation  for  this  problem  provided  by  Nakamura 
et  al.  (1985)  using  a  nodal  relaxation  scheme.  The  loading  and  geometry  for  this  example  are 
displayed  in  Fig.  5.  The  loading  is  applied  “suddenly”  at  the  top  and  bottom  edge  of  the  plate.  The 
crack  is  allowed  to  propagate  with  a  prescribed  motion.  In  particular,  the  crack  is  kept  stationary 
for  a  prescribed  amount  of  time  after  which  the  crack  is  made  to  advance  in  with  a  constant 
speed.  Because  of  this,  the  problem  solution  is  expected  to  reflect  the  discontinuous  nature  of 
the  crack  propagation  process.  Clearly,  the  problem  is  somewhat  artificial  in  that  the  crack  tip 
motion  is  a  given  of  the  problem.  However,  the  problem  is  mathematically  well  posed  and  it  is 
interesting  to  demonstrate  the  effectiveness  of  the  proposed  DGFEM  in  dealing  with  discontinuous 
loading.  For  the  time  interval  considered,  in  which  the  crack  tip  is  not  hit  by  reflected  waves,  the 
problem  affords  a  closed-form  solution  provided  by  Freund  (1990).  The  nondimensional  energy 
release  rates  calculated  by  the  formulation  developed  by  the  PI,  by  the  semidiscrete  method  used 
by  Nakamura  et  al.  (1985),  as  well  as  by  the  exact  solution  are  shown  in  Fig.  6  for  vc  =  0.4cs  and 
Fig.  7  for  vc  =  0.6cs,  vc  being  elastic  shear  wave  speed.  It  can  be  observed  that  the  numerical 
results  calculated  by  the  DGFEM  formulation  is  not  only  closer  to  the  closed-form  solution,  but, 
most  importantly,  it  captures  the  energy  release  rate  jump  discontinuity  caused  by  the  onset  of  the 
crack  motion. 

2.2.  DGFEM:  most  recent  developments 

The  result  discussed  above  were  obtained  with  a  FORTRAN  code  developed  from  scratch.  The 
code  in  question  was  rather  crude  and  not  easily  expandable  to  include  a  good  automatic  mesh 
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Figure  5:  The  geometry,  crack  position  and  boundary  conditions  of  for  the  results  in  Figs.6  and  7. 


Figure  6:  Energy  release  rate  estimates  for  vc  —  0Acs  and  step  size  At  =  l.OL/ci. 


refinement  capability.  Furthermore,  the  PI  intended  to  avoid  the  non-trivial  tack  of  creating  an 
efficient  program  with  error  estimation  and  adaptivity.  Therefore,  the  PI,  along  with  Scott  T.  Miller, 
a  Master’s  student  of  his,  decided  to  undergo  a  complete  restructuring  of  the  numerical  capability 
discussed  earlier,  so  as  to  rely  on  a  very  well  established  FEM  library.  After  some  searching,  it  was 
decided  that  any  new  FEM  development  was  to  be  based  on  the  “deal.ii”  FEM  library  developed 
and  maintained  by  W.  Bangerth,  G.  Kanschat,  and  R.  Hartmann  Bangerth,  2002;  Hartmann  and 
Houston,  2002a, b.*  deal.II  is  a  C-f~+  program  library  specifically  designed  to  support  adaptive 
finite  elements  and  error  estimation.  It  uses  state-of-the-art  programming  techniques  of  the  C-f-b 
programming  language  so  as  to  allow  one  to  use  complex  data  structures  and  algorithms  required 
for  adaptivity  from  ID  through  4D.  As  it  can  be  learned  from  the  supporting  literature  on  deal.ii, 

deal.II  emerged  from  work  at  the  Numerical  Methods  Group  of  the  University  of  Heidel- 
*For  detailed  information  concerning  the  deal.ii  library  see  <www.dealii.org>. 
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Figure  7:  Energy  release  rate  estimates  for  vc  =  0.6  cs  and  step  size  At  =  1.0L/ci. 


berg,  Germany,  which  is  at  the  forefront  of  adaptive  finite  element  methods  and  error 
estimators.  It  is  presently  developed  in  Heidelberg,  at  the  Institute  of  Aerodynamics 
and  Flow  Technology  of  the  German  Aerospace  Center  (DLR)  in  Braunschweig,  and  at 
the  Institute  for  Computational  Engineering  and  Sciences  (ICES,  formerly  TICAM)  in 
Austin,  Texas.  It  is  used  in  several  large  research  projects  around  the  world.  deal.II 
presently  has  about  220,000  lines  of  code  and  has  several  hundred  users  around  the 
world. 

Among  other  features,  the  deal.ii  library  is  extremely  thoroughly  documented  and  it  provides  fast 
algorithms  that  enable  you  to  solve  problems  with  up  to  several  millions  of  degrees  of  freedom, 
a  complete  stand-alone  linear  algebra  library  (including  sparse  matrices,  vectors,  Krylov  subspace 
solvers,  support  for  blocked  systems,  and  interface  to  other  packages  such  as  PETSc  and  METIS. 
Finally,  deal.ii  is  freely  available  under  the  Q  Public  Licence. 

The  most  recent  developments  of  the  PI  activities  consist  of  the  recasting  of  the  DGFEM  dis¬ 
cussed  earlier  using  deal.ii.  The  reformulation  in  question  has  been  successful  and  automatic  self 
refinement  is  now  part  of  the  numerical  capability  created  by  the  PL  In  the  process  of  recasting 
the  proposed  DGFEM  to  include  adaptivity,  it  was  discovered  that  the  formulation  used  until  then 
suffered  from  loss  of  unconditional  stability  when  using  completely  unstructured  space-time  grids. 
Hence,  the  recent  efforts  included  a  correction  of  the  formulation  and  the  derivation  of  a  proof 
of  unconditional  stability  of  the  new  formulation.  The  proof  in  question  has  been  presented  in  a 
Costanzo  and  Huang  (2004).  This  paper  includes  the  discussion  of  a  variant  of  the  proposed  DGFEM 
formulation  which,  in  the  absence  of  physically-based  dissipation,  is  energy  conserving  (although 
the  computer  implementation  of  this  variant  has  not  been  pursued  yet).  The  unconditional  stability 
of  the  new  formulation  has  been  also  numerically  verified.  Figure  8  presents  the  solution  in  terms 
of  a  simple  ID  problem  in  which  a  homogeneous  stretched  bar  initially  at  rest  is  suddenly  released. 
The  solution  was  obtained  a  completely  unstructured  space-time  grid.  The  total  energy  of  the  true 
solution  to  the  problem  is  conserved.  The  plot  of  the  total  energy  vs.  time  of  the  approximate 
solution  is  displayed  in  Fig.  9,  which  shows  a  (very  slight)  monotonic  decreasing  behavior  for  the 
energy,  thus  providing  a  numerical  verification  for  the  unconditional  stability  claim. 

Figure  10  presents  the  solution  for  the  same  phase  transition  problem  discussed  in  relation  to 
Fig.  2,  this  time  achieved  with  the  new  formulation  and  automatic  mesh  refinement.  Despite  the 
fact  that  half  as  many  elements  were  used  in  each  space-time  slab  (with  respect  to  the  calculation 
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Figure  8:  Plots  of  the  (nondimensional)  axial  internal  force  (stress)  in  a  homogeneous  bar  (shown 
in  the  inset  in  the  upper  left)  which  was  initially  at  rest  and  stretched.  At  time  t  =  0,  the  bar 
is  released  thus  causing  a  stress  pulse  to  propagate  to  the  left.  The  solid  black  line  denotes  the 
DGFEM  solution.  The  dashed  black  line  denotes  exact  solution,  whereas  the  gray  line  denotes 
the  solution  obtained  using  FEMLAB. 


Figure  9:  Plot  of  the  total  energy  normalized  with  respect  to  the  initial  total  energy  for  the  case 
described  in  Fig.  8.  The  solid  black  line  denotes  the  DGFEM  solution.  The  dashed  black  line 
denotes  exact  solution.  The  inset  plot  is  simply  is  a  magnification  of  the  main  plot  by  a  factor  of 
100  to  show  the  monotonic  decay  of  the  energy. 
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Figure  10:  The  dashed  line  represents  the  exact  solution.  The  solid  back  line  represents  the 
distribution  of  the  axial  internal  force  (stress)  at  various  time  instants  for  the  case  of  the  interface 
traveling  at  a  constant  speed  v,  with  v/c2  =  0.1.  C\jc2  —  0.5,  ki/k2  =  0.25,  and  n*(t  =  0)  =  0.1. 
Every  plot  has  the  same  horizontal  and  vertical  axes.  The  horizontal  axis  represents  the  non- 
dimensional  position  x *  along  the  bar,  while  the  vertical  axis  represents  the  non-dimensional  axial 
force  n*.  Each  space- time  slab  was  discretized  using  100  bi-quadratic  isoparametric  space-time 
elements  and  then  automatically  refined.  The  final  number  of  elements  per  slab  is  approximately 
200. 


in  Fig.  2),  the  solution  accuracy  remains  very  satisfactory.  Figure  11  displays  the  FEM  grid  that 
was  automatically  generated  over  the  space-time  solution  domain  using  automatic  refinement.  In 
this  figure,  the  automatic  refinement  captures  the  propagation  of  shock  wave  fronts  as  well  as  the 
propagation  of  the  interface.  Figure  12  shows  the  total  energy  vs.  time  plot  for  the  propagating 
interface  problem.  Note  that  the  dissipation  due  to  the  formulation  is  essentially  negligible  when 
compared  to  the  physical  dissipation  due  to  the  advancement  of  a  moving  boundary.  Finally,  the 
computational  capability  created  thus  far  was  used  to  predict  the  motion  of  an  interface  in  a  ID 
dynamic  solid-solid  phase  transition  process  using  as  evolution  law  of  the  type 


V(f  -  fcR ), 

0, 


for  /  >  fCR 
for  /  <  fcR, 


(2) 


where  is  the  velocity  of  the  interface,  r\  is  a  material  parameter,  /  is  the  energy  release  rate 
and  fcR  is  a  threshold  value  of  the  energy  release  rate  below  which  interface  propagation  does  not 
occur.  Formally,  the  structure  of  this  ID  problem  is  vary  similar  to  that  of  the  dynamic  propagation 
of  Griffith  crack  in  a  linear  elastic  material.  Figure  13  depicts  the  interface  trajectories  for  several 
values  of  the  r)  parameter  and  under  the  same  loading  conditions  of  the  problem  discussed  in  Fig.  10. 

The  results  reported  here  are  being  compiled  into  a  paper,  which  is  currently  in  preparation. 
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Figure  11:  Space-time  grid  automatically  generated  by  the  adaptive  refinement  for  the  solution  in 
Fig.  10. 


Figure  12:  Plot  of  the  nondimensional  total  energy  vs.  nondimensional  time  for  the  solution  in 
Fig.  10.  Both  exact  and  approximate  solutions  are  presented  in  the  plot. 

2.3.  Additional  developments 

Thorough  a  cooperation  with  Prof.  Jay  R.  Walton,  the  PI  is  working  on  a  formulation  of  dynamic 
fracture  in  unbounded  media  based  on  the  use  of  complex  variables  and  integral  transforms.  The 
PI  work  in  this  field  has  been  discussed  in  a  number  of  papers  preceding  the  beginning  the  present 
research  project  Costanzo  and  Walton,  1998a, b,  2002.  With  the  present  research  project,  the  PI 
has  begun  to  extend  the  solution  methods  already  developed  to  the  solution  of  a  Mode  III  transient 
dynamic  propagation  of  a  CZ  in  a  fully  coupled  thermo-elastic  solid.  The  progress  achieved  on 
this  component  of  the  project  consists  of  extending  the  computational  capability  already  created 
by  the  PI  to  solve  transient  (i.e.,  not  steady  state)  dynamic  fracture  mode  III  problems  with  a 
CZ  undergoing  fragmentation.  The  new  components  of  the  computational  capability  in  question 
regard  the  addition  of  adaptivity  as  well  as  being  able  to  deal  with  a  a  fragmented  CZ,  that  is, 
a  multiply  connected  domain.  The  complication  introduced  by  the  presence  of  fragmentation  is 
rather  substantial  from  a  computational  geometry  viewpoint.  The  progress  to  date  on  this  part  of 
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Figure  13:  Interface  trajectory  as  a  function  of  the  non-dimensional  time  t*.  From  left  to  right,  The 
various  curves  correspond  to  the  following  value  of  the  nondimensionalized  parameter  rj  in  Eq.  (2): 
20,  25,  30,  35,  40,  45,  50,  60,  75. 

the  project,  consists  in  the  creation  of  a  triangulation  capability  which.  can  in  indeed  cope  with  an 
evolving  and  fragmenting  CZ.  In  addition,  an  algorithm  for  the  use  of  product  integration  scheme 
to  deal  with  the  double- singular-integral  Neumann- to-Dirichlet  map  which  is  at  the  core  of  the 
numerical  scheme  in  question  has  also  been  formulated.  This  integration  scheme  guarantees  that 
all  of  the  numerical  integrations  to  be  carried  out  are  exact  for  the  interpolation  functions  selected 
to  represent  the  solution. 

Finally,  an  additional  development  is  the  extension  of  the  DGFEM  formulation  discussed  earlier 
to  the  solution  of  a  fully  coupled  thermo-elasto-dynamic  moving  boundary  problems.  This  compo¬ 
nent  of  the  project  has  begun  in  the  last  month  and,  while  a  few  results  have  been  already  obtained, 
a  more  careful  validation  of  the  numerics  is  needed. 

2.4.  Publications 

The  work  done  since  the  beginning  of  the  project  has  resulted  in  the  following  publications,  which 
acknowledge  AFOSR  funding: 

1.  F.  Costanzo  AND  H.  Huang,  “Proof  of  Unconditional  Stability  for  a  Single-field  Discontin¬ 
uous  Galerkin  Finite  Element  Formulation  for  Linear  Elasto-dynamics,”  Computer  Methods 
in  Applied  Mechanics  and  Engineering ,  to  appear,  2004. 

2.  H.  Huang  and  F.  Costanzo,  “On  the  Use  of  Space-time  Finite  Elements  in  the  Solution 
of  Elasto-dynamic  Fracture  Problems,”  International  Journal  of  Fracture,  to  appear,  2004. 

3.  Scott  T.  Miller,  “A  Space- Time  Finite  Element  Method  for  Second  Order  Hyperbolic 
Problems,”  Master  of  Science  thesis,  The  Pennsylvania  State  University,  August  2004. 

Mr.  Scott  T.  Miller  has  defended  his  thesis  work,  carried  out  under  the  supervision  of  the  PI,  in 
June  2004.  Mr.  Miller  is  going  to  continue  his  graduate  studies  in  the  Department  of  Theoretical 
and  Applied  Mechanics  at  the  University  of  Illinois — Urbana. 


3.  Interactions 

2002  Along  with  his  collaborator  Jay  R.  Walton  (Department  of  Mathematics,  Texas  A&M  Univer¬ 
sity),  the  PI  has  had  the  opportunity  to  give  a  seminar  to  the  GERC  in  April  2002,  illustrating 
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a  number  of  preliminary  results  obtained  on  the  use  of  the  DG  FEM  in  dynamic  fracture  of  lin¬ 
ear  elastic  isothermal  materials.  The  PI  intends  to  expand  his  collaboration  with  Dr.  Walton 
aimed  at  the  derivation  of  the  Neumann-to-Dirichlet  maps  for  coupled  linear  thermo-elastic 
dynamic  fracture  problems.  Furthermore,  while  at  GERC  in  April,  the  PI  had  the  oppor¬ 
tunity  to  meet  Dr.  Molly  Hughes  (Damage  Mechanisms  Branch  of  the  Ordinance  Division 
within  the  Munition  Directorate  at  AFRL/Eglin  AFB).  From  the  brief  discussion  with  her,  it 
appears  that  Dr.  Hughes  is  interested  in  the  Pi’s  project  and  possible  future  collaborations. 

2003  The  PI  has  been  contacted  by  Prof.  Barna  Szabo  (Professor  of  Mechanics,  Washington  Uni¬ 
versity)  for  exploring  the  full-scale  development  of  the  DG  FEM  scheme  proposed  herein. 
Furthermore,  the  PI  will  be  spending  two  weeks  visiting  with  Prof.  Walton  (October  16- 
November  1,  2003)  to  further  their  development  on  the  derivation  of  Neumann-to-Dirichlet 
maps  for  transient  thermo-elastic  fracture  problems. 

2004  The  PI  continues  to  interact  with  prof,  Jay  R.  Walton  on  matters  that  very  closely  related 
to  this  project.  Furthermore,  the  PI  is  still  in  coontact  with  Prof.  Szabo  for  exploring  the 
full-scale  development  of  the  DG  FEM  scheme  proposed  herein.  Finally,  the  PI  has  recently 
hired  as  a  post-doctoral  researcher  on  this  project  Dr.  Dinara  Khalmanova  who  received  her 
Ph.D.  in  mathematics  at  Texas  A&M  University  under  Prof.  Walton’s  supervision. 
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